The dynamics of spiral wave in the reaction-diffusion systems^[1–4] has been a hot topic of various scientific studies in various systems such as physical systems,^[5,6] chemical systems,^[7,8] myocardial tissues,^[9,10] the self-organization of slime mould,^[11] or neocortex.^[12] Experiment has evidenced that^[9] ventricular tachycardia could be the results of electrical spiral wave. The transition from ventricular tachycardia to ventricular fibrillation corresponds to the breakup of the spiral wave, in which an initiated single spiral wave breaks up into turbulence after several rotations. So it is of practical significance to study the behavior and control of spiral waves. Because of the potential applications, controlling spiral waves in the reaction-diffusion system has attracted a great deal of attention from scientists.^[13–17]

Currently, most studies of the dynamics of spiral wave in the heart only just consider the coupling among cardiomyocytes. It is well known that cardiac muscle mainly consists of myocytes and fibroblasts.^[18–20] Fibroblast content increases with its normal development and aging. Fibroblast proliferation and tissue remodeling are features of an aging heart. Myocytes are electrically coupled to other myocytes and fibroblasts by gap junction channels, forming a complex network of electrical impulse propagation. Although the propagation of electrical impulse is mainly accomplished by cardiac myocytes, fibroblasts can also play an active role in tissue electrophysiology. The experimental results^[21] show that fibroblasts of cardiac origin can synchronize the electrical activity between distant myocytes even if the length of fibroblast insert reaches . However, synchronization is accompanied with extremely large local conduction delays. Increasing insert length results in propagation failure due to the increase of activation delay. On the other hand, fibroblasts coupled to myocytes can act as a current sink in some cases. By imposing an electrical load, they can induce the regional shortening of action potential duration and affect myocardial excitability, and contribute to slow conduction.^[19] In ventricular tissue, fibroblast-based current sinks could cause a unidirectional block of the conduction.^[20] In addition, Nguyen et al.^[22] investigated whether myofibroblast-myocyte coupling can promote arrhythmia triggers by using the dynamic voltage clamp technique. They found that the coupling of myocytes to myofibroblasts could promote early after depolarizations (EAD) formation. Thus the fibroblast-myocyte coupling (FMC) mainly plays a role in inhibiting electrical impulse propagation. Hence, to improve our understanding of cardiac function in normal and pathological myocardium, it is important and necessary to extend homogeneous cell coupling to heterogeneous cell coupling.

In order to understand the effect of FMC on the cardiac electrical property, we propose a two-dimensional excitable system with an inhibitory coupling network. We investigate the effect of the inhibitory coupling on the dynamics of spiral wave by using the Bär–Eiswirth model.^[23] The numerical simulation results show that the inhibitory coupling has a great influence on the dynamics of spiral wave if the inhibitory coupling strength is big enough. Here we show that the inhibitory coupling can change the excitability of a system, the wave amplitude and duration of the excited phase, and even leads to the complete loss of system excitability and generates local abnormal oscillation in the system. Some new phenomena are observed in those systems with different parameters.

The outline of this paper is organized as follows. In Section 2 we introduce the mathematical model. In Section 3, we numerically investigate the effect of inhibitory coupling on the dynamics of spiral waves. The results of numerical simulations are presented. The results are summarized and discussed in the final section.

We use the Bär–Eiswirth model^[23] of excitable medium with the inhibitory coupling to explore the role of inhibitory coupling. The sites of the inhibitory coupling are uniformly distributed in a two-dimensional medium. The topology of the corresponding network (called the inhibitory network) will be discussed below. The fast variable u and slow variable v in the model obey the following differential equations:

All parameters and constants are dimensionless. The constant D is the self-diffusion coefficient, and the constants a, b, and ε are the system parameters. The last term on the right-hand side of Eq. (1a) is the inhibitory coupling term and is used to characterize FMC, d is the inhibitory coupling strength, and γ denote the position vector and the number of neighbors of the node in the inhibitory network, respectively, is the average value of the γ neighbor of u. It is obvious that the target state is compressed γ times compared with the self-diffusion term. Thus represents the inhibitory coupling. In our study, the positive parameters D, a, and b are fixed at D = 1.0, a = 0.84, and b = 0.07 to take advantage of the many known results. The results of Bär and Eiswirth^[23] show that the system governed by Eqs. (1a) and (1b) with d = 0 is excitable. The excitability of the system may be defined as the inverse of ε. Increasing ε can reduce the excitability of the system. Suitable initial conditions lead to rigidly rotating spiral waves (i.e., stable spiral waves) for . The spiral wave undergoes a transition from steady rotation to meandering when . Spiral waves break up into turbulence (i.e. spatiotemporal chaos) for .

In numerical simulation, the system size is taken as . The two-dimensional square medium is discretized into grid points. No-flux boundary conditions are used in all simulations. Equations (1a) and (1b) are integrated by the forward Euler method together with a second-order accurate finite-difference method with fixed time step and spatial step . We distribute an inhibitory coupling site for per m sites on the grid in the x and y directions respectively to generate an inhibitory network. Thus the system is discretized into a duplex network as shown in Fig. 1. The couplings in the horizontal and vertical direction are self-diffusive, whereas the couplings in the diagonal direction are inhibitory couplings. In the next section we investigate how the inhibitory coupling strength d and the distance between grid points linked by inhibitory coupling (i.e., m) affect the dynamics of a spiral wave in the duplex networks. Throughout the paper, spiral waves are generated by truncating planar waves.

Fig. 1.

	Figure Option ViewDownloadNew Window
	Fig. 1. (color online) Schematic diagram of the duplex networks. Black and red links represent self-diffusive and inhibitory couplings, respectively. The red grids correspond to the inhibitory network with m = 2.

Since the behavior of spiral waves depends sensitively on parameters m, d, and ε, we take those parameters as adjustable parameters of the model. We let , and vary d in a broad range so that the systems can fully exhibit their dynamical behaviors. In the following we will study systematically the dynamics of spiral wave on the duplex networks.

3.1. Inhibitory network with shortest range connection

We next show how the inhibitory coupling affects the dynamics of spiral waves in duplex networks with m = 1, which corresponds to the pathological myocardium. We find that the inhibitory coupling directly affect the excited phase duration (EPD), wavelength λ of spiral wave, wave amplitude (i.e., maximal value of variable u) and excitability of the system. Increasing inhibitory coupling strength can induce the phase transition of the system.

When ε = 0.01, the wavelength λ of spiral wave will increase slowly as d is increased up to d = 3.3, at which point the spiral wave starts to meander. After the transition point, the wavelength increases quickly as d is further increased. The variation of wavelength with d is shown in Fig. 2. It is observed that the wavelength increases with inhibitory coupling strength increasing due to the decreased excitability of the system.

Fig. 2.

	Figure Option ViewDownloadNew Window
	Fig. 2. Patterns of spiral waves for d = 1 (a), 3.5 (b), and 3.8 (c), with parameters m = 1 and ε = 0.01. The gray level is proportional to u variable.

The inhibitory coupling can result in the decreases of wave amplitude and EPD. Comparing with the scenario without inhibitory coupling, the EPDs at u = 0.1 (EPD₉₀) for the parameters of Figs. 2(a)–2(c) decrease by 15.5%, 24.3%, and 19.4%, respectively. Figures 3(a) and 3(b) show the dependence of and T on d for and m = 1, respectively. From Fig. 3(a) we can see that when d = 0, the wave amplitude is equal to 1. When the inhibitory coupling strength d is increased, will decrease monotonically until decreases to 0.727 at , beyond which the inhibitory coupling completely suppresses the formation of a spiral wave. A self-sustaining planar wave is generated by truncating a planar wave (see Fig. 4). From Fig. 4, it is observed that the planar wave moves along the border of the system. After the planar wave turns a corner, it will rotate around the free end, which leads to the decrease of its length. Its length is restored to its original length after rotating. When d further increases, the maximum length of planar wave will decrease until the planar wave shrinks into a spot wave. The corresponding reaches its minimum . When , there is no wave propagation due to conduction failure. We refer to and as the critical inhibitory coupling strength and the minimal wave amplitude, respectively. That is the reason why is equal to 0 for . Figure 3(b) shows that the period of spiral wave increases slowly before d reaches 3.3, after which point the period increases rapidly until the system is not able to propagate a wave. The above results reveal that the increasing inhibitory coupling strength can gradually reduce or even completely suppress the excitability of the excitable system. The maximal decrease in wave amplitude can reach 28.1%. Maleckar et al. obtained similar results by using an electrophysiological model.^[24] They found that when a human atrial myocyte couples to two active 1 fibroblasts for low gap junctional conductance, FMC results in the dramatic decrease of myocyte excitability. The action potential amplitude and duration (APD₉₀) of the atrial myocyte can reduce 12.4% and 17.3% respectively as compared with those without FMC. In contrast, coupling to active 2 fibroblasts causes myocyte APD₉₀ to increase by 16.3%. The action potential amplitude reduces 34% for high gap junctional conductance. The electric current between fibroblasts and a myocyte can lead to complete loss of myocyte excitability.

Fig. 3.

	Figure Option ViewDownloadNew Window
	Fig. 3. Plots of (a) wave amplitude at a grid point and (b) the period of spiral wave T versus the inhibitory coupling strength d for m = 1 and ε = 0.01.

Fig. 4.

	Figure Option ViewDownloadNew Window
	Fig. 4. Movement of planar wave at m = 1, ε = 0.01, d = 4.05, and t = 2300 (a), 2340 (b), 2380 (c), and 2420 (d).

When ε is increased, the numerical results show that the critical inhibitory coupling strength monotonically decreases. For example, when ε increases to 0.03, decrease to 1.2. When d is increased, the amplitude of spiral wave monotonically decreases for a given value of ε. The transition from stable spiral wave to meandering spiral wave, the transition from meandering spiral wave to planar wave transition and the transition from planar wave to no wave successively appear for ε = 0.02 and ε = 0.03. For ε = 0.03, we find that EPD decreases at first. Then EPD increases again with d increasing further. The maximal increment of EPD₉₀ reaches 13.4%.

When ε = 0.035, the spiral wave exhibits the unique meandering as illustrated in Fig. 5, which shows that the spiral wave first breaks up in the spiral center after about three rotations, leading to the tip of the spiral wave disappearing. Then the spiral curve starts to drift until the free end creates a new spiral. The process repeatedly appears, leading to the meandering of the spiral wave. In order to understand the mechanism behind the spiral wave breakup, we plot the time evolution of u as shown in Fig. 6. It is observed that u exhibits abnormal oscillations at the excited phase of u, which is generally referred to as early after depolarizations (EADs). It is obvious that the local abnormal oscillations of the system result in spiral wave breakup.

Fig. 5.

	Figure Option ViewDownloadNew Window
	Fig. 5. Time evolutions of meandering spiral wave at m = 1, ε = 0.035, d = 0.9, and t = 184 (a), 198.4 (b), 200.8 (c), 206.8 (d), 213.2 (e), and 229.6 (f).

Fig. 6.

	Figure Option ViewDownloadNew Window
	Fig. 6. Time evolution of u at a grid point inside the spiral center shown in Fig. 5(b) for m = 1, ε = 0.035, and d = 0.9.

For ε = 0.04 and ε = 0.05, the transition from meandering spiral wave to planar wave vanishes. Spiral wave is stable when dis small. As d is increased, turbulence and a no-wave state arise successively in the system. Turbulence mainly comes from two different ways: one comes from the Doppler effect (i.e., a spiral wave develops into turbulence first near the spiral core as illustrated in Fig. 7); the other comes from the backfiring instability. The backfiring is the re-excitation behind the pulse of activation, and it generates the emission of a pulse in the opposite site direction as illustrated in Fig. 8. In Fig. 8, the ring and the short lines contacting the longest line are generated by the backfiring effect, which is generated by the abnormal oscillation of u. The two types of turbulences are generally called turbulences 1 and 2 (T1 and T2), respectively.

Fig. 7.

	Figure Option ViewDownloadNew Window
	Fig. 7. Development of turbulence at m = 1, ε = 0.04, d = 0.6, and t = 50 (a), 60 (b), 80 (c), and 500 (d).

Fig. 8.

	Figure Option ViewDownloadNew Window
	Fig. 8. Time evolutions of turbulence at m = 1, ε = 0.04, d = 0.9, and t = 45 (a), 80 (b), 120 (c), and 500 (d).

When , spiral waves in the system without inhibitory coupling meander regularly. When d is slightly increased, the meandering spiral wave will break up into turbulence. When , the system is in a no-wave state. For example, for ε = 0.065, meandering spiral waves can be seen only at , whereas for ε = 0.069 spiral waves break up into turbulence even if a smaller d is applied. For , the increase of d only results in the transition from turbulence to the no-wave state.

3.2. Inhibitory network with long range connection

Now we investigate the dynamics of spiral waves in the duplex networks with . The numerical results show that the increase of inhibitory coupling strength d also leads to the transition to meandering spiral wave, turbulence or no-wave for a given ε and m. However, the coupling strength d relating to the transition increases with m increasing. When , the transition to no-wave vanishes. For , the inhibitory coupling has hardly any influence on the dynamics of the system. These results suggest that the effects of FMC on myocyte and the excitability of the system will become more and more weak as the fibroblast content decreases.

When the parameters are chosen at the boundary between turbulent regions T1 and T2, we observe the Eckhaus instability of spiral wave for as shown in Fig. 9. The wide circular wave shown in Fig. 9 is generated by the local abnormal oscillations of the system. The abnormal oscillation-induced backfire in turn results in an inward spiral wave. After the circular wave moves out of the system boundary, the spiral wave begins to break up due to the Eckhaus instability, and it develops into turbulence finally.

Fig. 9.

	Figure Option ViewDownloadNew Window
	Fig. 9. Development of turbulence at m = 3, ε = 0.07, d = 3, and t = 50 (a), 70 (b), 80 (c), 90 (d), 100 (e), and 500 (f).

In order to characterize the dynamical behavior of the system with the inhibitory coupling, a phase diagram in the d–ε plane is drawn in Fig. 10. We can see clearly from Fig. 10 that: (i) upon increase of ε, decreases gradually for a given m. The area of the no-wave region decreases gradually as m is increased until it vanishes at a larger m. For , the phenomenon of the no-wave state does not occur. (ii) There are two separate regions of the meandering spiral wave for . However, in the case of m = 1, the two meandering regions are small, thus the meandering region relating to the large d is not shown in Fig. 10(a). (iii) There are two kinds of transitions to the no-wave state. They are the transitions from planar wave and turbulence to no wave, respectively. In addition, there are the transitions from stable and meandering spiral wave to turbulence. Those transitions can be achieved by changing at given . (iv) When m is large enough, the system will experience different scenarios as ε is increased, the scenarios are similar to those of the system without inhibitory coupling. In other words, the inhibitory coupling has hardly cause any influence on the dynamics of a system when m is large enough.

Fig. 10.

	Figure Option ViewDownloadNew Window
	Fig. 10. Phase diagrams revealing regions of different wave forms: stable spiral waves (i.e., rigidly rotating spiral waves) (S), meandering spiral waves (M), turbulence 1 (T1, cf. Fig. 7), turbulence 2 (T2, cf. Fig. 8), no wave (N), for m = 1 (a), 2 (b), 3 (c), 4 (d), 5 (e), and 6 (f).

We investigated the dynamics of spiral waves in the duplex networks with the inhibitory coupling. We find that the complex effects of the inhibitory coupling on the dynamics of spiral wave depend on the specific properties of the system. Both EDP shortening and prolongation as well as wave amplitude decrease are observed. For the inhibitory network with m = 1, the maximal decrement (maximal increment) in EDP₉₀ reaches 24.3% (13.4%). The maximal decrement in wave amplitude approaches 28.1%. Upon increase of d the ability for a system to propagate a wave is lost, causing a stable spiral wave to develop into a meandering spiral wave or turbulence, even leading to conduction failure. When parameters are appropriately chosen, the inhibitory coupling can induce the abnormal oscillation of u. The abnormal oscillation can either induce the unique meandering of a spiral wave, or cause the transition to a turbulence state. The results are essentially consistent with those of the electrophysiological model and experiment. In addition to those results, we have some important findings. We find that the inhibitory coupling can generate a self-sustaining planar wave and induce the Eckhaus instability of a spiral wave in the excitable system with the inhibitory coupling, although an Eckhaus instability of a spiral wave generally occurs in an oscillatory system. The transitions from meandering spiral wave to planar wave and to inward spiral wave are observed. The inward spiral wave is generated by the local abnormal oscillation of the system. These results will help us to understand the phenomena in the heart.